Poniżej przedstawiam kod MATLABa obliczający problem adwekcji-dyfuzji metodą elementów skończonych ze stabilizacją metodą Streamline Upwind Petrov-Galerkin (SUPG), dla modelowego problemu Erikksona-Johnsona.
Przypomnijmy najpierw problem adwekcji-dyfuzji, problem modelowy Erikksona-Johnsona (opisany na przykład w dokumencie [1]), zdefiniowany na obszarze kwadratowym \( \Omega = [0,1]^2 \) w następujący sposób: Szukamy funkcji \( u \) takiej że
\( a(u,v)=l(v) \forall v \) gdzie
\( a(u,v) =\int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y) }{\partial x } dxdy + \int_{\Omega} \beta_y(x,y) \frac{\partial u(x,y)}{\partial y } dxdy \\ +\epsilon \int_{\Omega} \frac{\partial u(x,y) }{\partial x} \frac{\partial v(x,y)}{\partial x } dxdy +\epsilon \int_{\Omega} \frac{\partial u(x,y)}{\partial y } \frac{\partial v(x,y)}{\partial y } dxdy \)
\( l(v) = \int_{\partial \Omega } f(x,y) v dxdy \)
\( f(x,y)=sin(\pi y)(1-x) \) jest rozszerzeniem warunku brzegowego Dirichleta na cały obszar, natomiast
\( \beta = (1,0) \) reprezentuje wiatr wiejący z lewej strony na prawą, natomiast \( \epsilon = 10^{-2} \) oznacza współczynnik dyfuzji. Wielkość \( Pe=1/ \epsilon = 100 \) nazywa sie liczbą Pekleta, i definiuje ona wrażliwość problemu adwekcji-dyfuzji.
W linii 848 podawany jest współczynnik dyfuzji \( \epsilon = 10 ^{-2} \). Jego odwrotnośc to tak zwana wartość Pekleta (Peclet number)
\( Pe=\frac{1}{\epsilon} \).
W liniach 849 i 850 zaszyte są stałe metody SUPG.
Następnie tworzony jest wektor węzłów (zawsze liczby naturalne liczone od zera) \( knot\_x = [0 \quad 0 \quad 0 \quad 1 \quad 2 \quad 2 \quad 2] \) dla przestrzeni aproksymacyjnej,
oraz wektor odpowiadających im punktów wzdłuź osi x, z przedziału od 0 do 1, \( points\_x = [0 \quad 0.5 \quad 1] \) w których umieszczana będzie baza funkcji B-spline służaca do aproksymacji rozwiązania problemu Erikksona-Johnsona. Wektory węzłów dla przestrzeni aproksymującej i punkty na których rozpinamy wektory węzłów, definiujemy osobno dla osi x oraz dla osi y.
Kod może zostać uruchomiony w darmowym środowisku Octave.
Kod uuchamia się otwierając go w Octave oraz wpisując komendę
\( advection\_supg\_adapt \)
Po chwili obliczeń kod otwiera dodatkowe okienko i rysuje w nim rozwiązanie numeryczne oraz dokładne.
Kod oblicza błąd reziduum
\( Residual norm: L2 0.015808, H1 0.100417 \)
iKod oblicza również normę L2 i H1 z rozwiązania
\( Error: L2 31.41 procent, H1 38.74 procent \)
i porównuje do normy L2 i H1 z rozwiazania dokładnego.
\( Best possible: L2 2.06 procent, H1 28.73 procent \)
Widać że w celu poprawienia dokładności rozwiązania, konieczne jest zagęszczanie siatki.
Zadanie 1: Jednorodne zwiększenie siatki
Treść zadania:
Proszę zwiększyć liczbę punktów do pięciu w kierunku x i y, zachowując stopnie wielomianów B-spline. Proszę sprawdzić jak operacja ta poprawia dokładność rozwiązania
podobnie w kierunku y.
Error: L2 19.64procent H1 45.47procent
Best possible: L2 1.38procent H1 51.91procent
Zadanie 2: Adaptacyjne zwiększenie siatki
Treść zadania:
Proszę zwiększyć adaptacyjnie liczbę punktów w kierunku prawego brzegu [0 0.5 0.9 0.95 1] tam gdzie znajduje się warstwa przybrzeżna, zachowując stopnie wielomianów B-spline. Proszę sprawdzić jak operacja ta poprawia dokładność rozwiązania
Zadanie 3: Adaptacyjne zwiększenie siatki trial i test poprzez łamanie na pół elementu najbliżej osobliwości dla dużej wartości Peclet number
Treść zadania:
Proszę zmodyfikować problem tak aby liczba Pecleta wynosiła milion. Proszę wystartować od siatki [0 0.5 1] w kierunku x oraz [0 0.25 0.5 0.75 1] w kierunku y, i zwiększać adaptacyjnie liczbę punktów w kierunku prawego brzegu tam gdzie znajduje się warstwa brzegowa, zachowując stopnie wielomianów B-spline w przestrzeni trial i test. Proszę sprawdzić jak operacja ta poprawia dokładność rozwiązania
Rozwiązanie:
Wprowadzamy \( \epsilon=10^{-6} \) oraz modyfikujemy siatkę w kierunku x zgodnie z sekwencją:
.